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Abstract 


A pore network model of the gas diffusion layer (GDL) in a polymer electrolyte membrane fuel cell is developed and validated. The model 
idealizes the GDL as a regular cubic network of pore bodies and pore throats following respective size distributions. Geometric parameters of the 
pore network model are calibrated with respect to porosimetry and gas permeability measurements for two common GDL materials and the model 
is subsequently used to compute the pore-scale distribution of water and gas under drainage conditions using an invasion percolation algorithm. 
From this information, the relative permeability of water and gas and the effective gas diffusivity are computed as functions of water saturation using 
resistor-network theory. Comparison of the model predictions with those obtained from constitutive relationships commonly used in current PEMFC 
models indicates that the latter may significantly overestimate the gas phase transport properties. Alternative relationships are suggested that better 
match the pore network model results. The pore network model is also used to calculate the limiting current in a PEMFC under operating conditions 
for which transport through the GDL dominates mass transfer resistance. The results suggest that a dry GDL does not limit the performance of a 


PEMFC, but it may become a significant source of concentration polarization as the GDL becomes increasingly saturated with water. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Polymer electrolyte membrane fuel cells (PEMFCs) are a 
promising energy conversion technology. However, there are 
still several technological difficulties that must be overcome 
before they can be commercialized. One of the main challenges 
is to achieve effective water management inside the cell, since 
the presence of water can be both detrimental and beneficial to 
PEMFC performance and durability. A highly humidified envi- 
ronment is preferred in the cell to maintain membrane hydration 
and conductivity. Excess humidity, however, results in conden- 
sation and blockage of pores in the electrode backing or gas 
diffusion layer (GDL). These effects are complicated by the fact 
that water is a product of the oxygen reduction reaction in the 
cathode compartment. At high current densities, the increased 
rate of water production can lead to liquid water formation and 
flooding of the GDL. An additional difficulty is that the environ- 
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mental conditions inside the cell, such as temperature, pressure 
and gas compositions, can vary widely over the active area of 
a cell [1,2]. As a result, ideal humidity conditions may exist in 
one location while liquid water may form elsewhere. Clearly, 
understanding of the formation, behavior and movement of liq- 
uid water inside the porous components of the PEMFC is of 
great importance. 

A large number of multiphase flow models have recently 
appeared in the literature that attempt to address the problem 
of liquid water behavior in the cathode and its impact on mass 
transfer in a PEMFC [3-10]. The models presented to date 
are exclusively based on continuum descriptions of flow and 
transport, which require knowledge of constitutive relationships. 
These include the dependences on water saturation of the relative 
permeability, effective diffusivity and air—water capillary pres- 
sure. At present, GDL-specific experimental data on gas or liquid 
phase relative permeability are scarce, the effective diffusivity 
has been estimated only from numerical models [11] and only 
recently have air—water capillary pressure data been made avail- 
able [12]. As a result, many of the necessary relationships and 
parameters incorporated in elaborate multiphase transport mod- 
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Nomenclature 

a exponent in Eqs. (18) and (20) 

A area of lattice normal to flow direction (m?) 

b conduit size (radius) (m) 

c concentration (mol m~?) 

df fiber diameter (m) 

D diffusion coefficient (m? s7!) 

ga diffusive conductivity (m3 s7!) 

gh hydraulic conductivity (m? Pa~! s7!) 

K permeability (m7) 

l length of pore network domain (m) 

L length of conduit (m) 

Le lattice constant (m) 

L throat length (m) 

n diffusion rate through a pore conduit (mol s7!) 

N diffusion rate through network (mol s7!) 

Pc capillary pressure (Pa) 

q flow rate through a pore conduit (m? s7!) 

Q flow rate through network (m3 s7!) 

s saturation 

V volume (m?) 

X lattice size in x-direction (in-plane) no. of pores 

Y lattice size in y-direction (in-plane) no. of pores 

Z lattice size in z-direction (through-plane) no. of 
pores 

Subscripts 

b bulk 

B species B 

CH gas channel 

CL catalyst layer 

eff effective 

in inlet 

nwp non-wetting phase 

max maximum 

out outlet 

p pore 

P phase 

r relative 

t throat 

T total 

wW water 

wp wetting phase 

x x-direction (through-plane) 

y y-direction (in-plane) 

Z z-direction (in-plane) 

Greek symbols 

a throat constriction factor 

B spatial correlation distance 

X random number in Weibull distribution 

ô GDL thickness (m) 

E porosity 

y surface tension (N m7!) 

n filling exponent 


K parameter in Weibull distribution 

À parameter in Weibull distribution 

u viscosity (Pa s) 

0 contact angle (radians) 

o conducting phase 

Superscripts 

g value at pore breakthrough pressure 
m exponent used in Eq. (14) 

n exponent used in Eq. (15) 


els remain uncertain and application of these models to different 
GDL materials is questionable. 

An alternative approach to modeling multiphase transport 
processes in GDL materials is pore network modeling. This 
approach has a long history in the study of porous media of geo- 
logic origin (soil and rock) [13—16]. The basis of this approach 
is a mapping of a complex pore space continuum onto a regular 
or irregular lattice of sites and bonds. To derive a geometrical 
model it is usually assumed that the pore space can be conceptu- 
ally partitioned into a collection of pore bodies communicating 
through local constrictions termed pore throats. Model pore net- 
works are thus constructed by assigning pore and throat sizes 
to the lattice sites and bonds, respectively. Simplifying assump- 
tions regarding the shape of pores and throats are invariably 
made to facilitate the computation of capillary and transport 
characteristics of the pore network elements [17]. Pore network 
models are ideally suited for the simulation of low-capillary 
number (quasi-static) immiscible displacement using percola- 
tion concepts [13]. A main advantage of pore network models 
is that they account explicitly for pore-level physics and pore 
space geometry/topology. Prediction of various macroscopic 
transport and capillary properties of porous media is relatively 
straightforward if the geometric, topological and correlation 
properties of the porous microstructure are properly specified. 
The task of extracting this information is, however, non-trivial, 
typically requiring extensive characterization of 3D volume 
data [18]. 

The present work outlines the development of a pore net- 
work model to study multiphase transport in GDLs. This is the 
first attempt to deploy pore network modeling for the study 
of the gas diffusion layer of a PEMFC, although Thompson 
[19] has applied a pore network modeling approach to conven- 
tional paper. Numerous modifications are made to the traditional 
pore network modeling framework in order to account for the 
unique geometric aspects of fibrous GDLs. In the absence of 
3D volume data for the GDL materials studied, the network 
parameters are obtained by calibration to experimental gas per- 
meability and drainage capillary pressure data. The model is 
then used to simulate multiphase transport scenarios of inter- 
est to PEMFC operation, such as the diffusion of gas through 
a partially water-filled GDL and the convective flow of gas and 
water under conditions of partial water saturation. Results are 
presented for two typical GDL materials for which the neces- 
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sary experimental information is available. Finally, calculations 
of limiting current densities are performed by placing typical 
fuel cell boundary conditions on the network model and calcu- 
lating the mass transfer flux through partially saturated GDLs to 
the catalyst layer. 


2. Model development 
2.1. Materials to be modeled 


In this work, the porous networks of two different types of 
carbon paper are modeled. Fig. 1 shows micrographs of SGL 
Sigracet® 10BA and Toray 090. Toray 090 has a mostly 2D 
structure with linear fibers arranged in layers in the plane of 
the paper. SGL 10BA has a more 3D structure with intertwined, 
curved fibers. Physical properties of each material are listed in 
Table 1. 


2.2. Pore network construction 


One of the distinguishing features of GDLs is that they pos- 
sess a very high porosity, which can range from 0.75 to above 
0.90, meaning that GDLs are predominantly void space. More- 
over, there is little constriction between pores, creating a highly 


Table 1 

Physical properties of GDL materials 

Property SGL 10BA Toray 090 
Thickness (8) (um) 390 290 

Total porosity (£) 0.88-0.90 0.78-0.80 
Fiber diameter (df) (wm) 9 9 
Permeability (Ky) (m?) 57 x 107!? 15 x 107!? 
Permeability (Ky) (m°) 45 x 10-2 15x107” 
Permeability (K;) (m?) 37 x 107! 9.0 x 107!? 


open structure. Fig. 2 shows a cross-sectional slice obtained from 
a simple solid model of a GDL. With such small solid phase frac- 
tion, it is difficult to define distinct pore bodies or to identify pore 
throats. This situation is quite different from the one encountered 
in rocks and soils, for which pore bodies and pore throats can 
be intuitively delineated in images of the pore space. 


2.2.1. Pore and throat size distributions 

The pore network model developed here for GDLs is based on 
the one described by Ioannidis and Chatzis [17] and Chang and 
Ioannidis [20]. The pores are modeled as nodes on a regular cubic 
lattice, interconnected with throats. The pores are idealized as 
cubic bodies and the throats are treated as ducts of square cross- 


Fig. 1. SEM micrographs of GDL materials modeled in present study. (a) Toray 090 and (b) SGL Sigracet 10BA. (i) 100x and (ii) 1000x magnification. The fiber 


alignment in the SGL10BA sample is apparent in (b, i). 
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(a) (b) 


Fig. 2. Slices of a simulated GDL solid model. (a) In-plane view, (b) through- 
plane view. Both views show 10 pm thick sections. The model was generated 
by placing fibers with a random location and in-plane rotation, then applying an 
out-of-plane rotation with angles normally distributed around 0° with a standard 
deviation = 1. 


section. This arrangement is shown in Fig. 3 with the relevant 
dimensions labeled. The use of square pores is convenient in 
order to achieve sufficiently high porosities and to qualitatively 
describe the presence of corners and crevices in the pore space. 
The pore network is constructed by assigning pore body sizes 
from a truncated Weibull cumulative distribution: 


bp,i = AL[— In (1 — XXmax)] 7! + bmin (1) 


where bp, is the radius of the ith pore, x is a random number 
between 0 and 1, Xmax (<1) scales the random number and trun- 
cates the upper end of the distribution to prevent excessively 
large pores from being generated, bmin is the minimum pore 
radius and A and «x are adjustable parameters that control the 
location/spread and shape of the distribution. A Weibull distribu- 
tion is used since it is highly versatile and mathematically simple 
[17], containing only two adjustable parameters. These features 


| ; y C] Pore Body [| Pore Throat Gy Solid 


Fig. 3. Schematic diagram of two neighboring pore bodies and connecting 
throat. Throat size (b;) is proportional to the size of the smaller of the two con- 
necting pores (b;=abp). Throat length (L+) is equal to the difference between 
the pore body sizes (bp) and the center-to-center distance between pores (Lc). 


are advantageous when the pore size distribution is adjusted to 
calibrate the model as described in Section 3.1. 

Once pore sizes are assigned, throat sizes are assigned by 
assuming that the size of each throat is equal to the size of 
the smallest of the two adjacent pores. This throat assignment 
scheme is chosen because it allows for minimum constriction 
between pore bodies, creating a highly open structure charac- 
teristic of GDLs. Fig. 4a shows the construction of the lattice 
with pores and throats identified. Fig. 4b shows only the void and 
solid space of the same lattice. The open nature of the pore space 
obtained by this method of throat size assignment is apparent. 

The length of each throat is calculated as the difference 
between the lattice constant Lc and the size of the two connecting 
pores. The lattice constant is the spacing between pore centers 
and is adjusted to match the porosity of the network model to 
the known porosity of the material. This is discussed further in 
Section 3.2. Consequences of this size assignment scheme are 
that throats and pores have similar size and their volume cannot 
be neglected in the calculation of the total lattice volume. In 
fact, a throat is actually an extension of the pore body to which 
it is attached and the lattice is basically an assembly of pores 
connected directly to pores. 

It should be clear that the aforementioned description is by 
no means an attempt to reproduce the actual geometry of GDL 
pore space. What is sought instead is to endow the pore network 
model with sufficient flexibility to reproduce experimental mea- 
surements of capillary pressure and gas permeability (in-plane 
and through-plane). Obviously, a better way to construct the 
pore network would be to extract its geometric and topologi- 
cal properties from experimental 3D volume data of the GDL 
materials. 


2.2.2. Spatial correlation of pores sizes 

One of the key features included in the model is spatial cor- 
relation of pore sizes. A highly porous material such as a GDL 
contains regions of extended continuous void space with no solid 
to mark distinct boundaries between pore bodies. In terms of the 
pore network model, these regions are analogous to multiple 
neighboring pores of similar size. Imposing spatial correlation 
of pore sizes in the model results in pores of similar size being 
placed next to each other in the lattice. These pores are invaded 
by the non-wetting phase at similar capillary pressures and offer 
similar resistance to fluid flow, therefore acting as a single, large 
pore. The effect of introducing spatial correlation of pores into 
the model is to increase the permeability of the network by 
more than 20% and bring it more in line with measured val- 
ues. Experience has shown that without spatial correlation, it is 
very difficult to match both the experimental permeability and 
the capillary pressure curves, since both are dependent on pore 
size distribution (see further discussion in Section 3). 

Spatial correlation also partially accounts for the observed 
directional anisotropy in the permeability tensor [21]. When 
pores are correlated in certain directions, the permeability along 
these directions is increased. It was found that correlating pores 
in the direction of fiber alignment helped to create the observed 
anisotropy trends. For instance, since the fibers of Toray 090 
are aligned in the x-y plane, correlation of neighboring pores in 
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Fig. 4. 2D schematic diagram of pore network construction. (a) Relationship between pores, throats and solid. (b) Structure in terms of void and solid space. 


this plane, but not in the through-plane (z-direction), produces 
the correct trend. This is summarized with the notation [B,, By, 
B.]=[1, 1, 0] where £ is the correlation distance. The fibers 
in SGL 10BA are also predominantly aligned in the x-y plane, 
but have additional directional alignment in the x-direction. The 
use of correlation distances [B,, By, Bz] =[2, 1, 0] partially repro- 
duces the observed anisotropy. Fig. 5a shows a structure obtained 
using a field of random, uncorrelated numbers, whereas Fig. 5b 
and c show the structures obtained when the correlations [1, 1, 
0] and [2, 1, 0], respectively, are imposed. 

Anisotropy can also be created in the model by constricting 
throat sizes along specific directions. In addition to the imposi- 
tion of spatial correlation, a small amount of throat constriction 
was necessary to completely match the experimentally observed 
anisotropy in permeability. Throats were uniformly constricted 


according to the expression: 
bij = æbp,i (2) 


where bij is the size of the throat connecting pores i and j, 
bp; is the size of pore i with bpij<bpj and œ is the throat 
constriction factor. The throat constriction factor is direction 
dependent and described with the notation [a,, ay, œz]. In gen- 
eral it was necessary to constrict throats slightly (S—10%) in the 
direction perpendicular to the axis of fiber alignment. For Toray 
090 throats were constricted in the through-plane z-direction 
according to [&x, ay, a] =[1, 1, 0.9]. In SGL 10BA, the fibers 
are aligned in the x-y plane with some additional alignment 
in the x-direction. Accordingly, throat constriction factors [a,, 
dy, az] =[1, 0.95, 0.95] were used. Constricting throats in this 


Fig. 5. Examples of spatially correlated random fields. (a) Uncorrelated field. (b) Correlated field used to model Toray 090 with correlation distances [1, 1, 0] in the 
x, y and z directions (z-direction not shown). (c) Correlated field used to model SGL 10BA with correlation distances [2, 1, 0] in the x, y and z directions (z-direction 


not shown). 
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way is consistent with the structure of GDLs since flow in the 
cross-fiber direction is more obstructed. 


2.3. Capillary pressure 


All pore throats and pore bodies in this model are assumed to 
be of square cross-section. The capillary pressure, Pc, required 
for a non-wetting fluid to penetrate a throat of square cross- 
section is estimated by the Young—Laplace equation: 


Po = —2y cos 6 (>) (3) 
bi 
where y is the surface tension, 6 is the contact angle and b; 
is the radius of the largest circle that can be inscribed in the 
square capillary. Contact angles in GDL materials are not easily 
determined. The contact angle on simple carbon is highly vari- 
able [22]. In a previous study [12], an experimental procedure 
was described for estimating the microscopic contact angle of 
mercury on GDL fibers by measuring the macroscopic contact 
angle of a sessile drop on the GDL surface and correcting for 
the porosity and roughness of the surface. It was estimated that 
the microscopic mercury contact angle could be as low as 110°. 
In the present work, an angle of 115° was used. We measured 
macroscopic contact angles of water on the two GDLs of this 
study and adopted the same procedure to estimate their corre- 
sponding microscopic contact angles. Table 2 lists the values so 
obtained. We note that GDL materials containing carbon and 
Teflon are expected to have non-uniform wettability, although 
no data are presently available to quantify this expectation. The 
pore network detailed here can be modified to accommodate 
non-uniform contact angles. 


2.4. Late pore filling 


In reality, pore geometry is more complex than any simple 
geometric shape, albeit angular, can describe. Unresolved length 
scales due to the presence of cracks, corners, crevices and inter- 
stitial regions at fiber—fiber contact points amount to pore space 
from which the wetting phase is displaced at capillary pressures 
higher than corresponding to first entry of the non-wetting phase 
into any pore in the network. To account for the gradual drainage 
of the wetting phase from such small scale features, we employ 
the following expression [20]: 


P* n 
Swp = stp( 75) , Pe> PŽ (4) 

Pe 
Table 2 
Fluid properties 
Fluid Surface tension (N m7!) Contact angle* (°) 

SGL 10BA Toray 090 

Mercury 0.480 115 115 
Water 0.072 100 98 
Octane 0.022 0 0 


è Measured through the liquid phase. 


where n is the filling exponent, swp is the wetting phase saturation 
of a given pore at capillary pressure Pc, and Sp is the wetting 
phase saturation of the same pore at the capillary pressure, P%, 
corresponding to first entry (breakthrough) of the non-wetting 
phase. The parameters n and Sp are adjustable. Late pore filling 
enables smaller scale features to affect the capillary pressure 
behavior of the network without explicitly including them as 
individual pores. This treatment was found to be necessary to 
correctly model the experimental capillary pressure curves. 


2.5. Drainage simulation 


The process considered by the present model is the drainage 
of a wetting phase by slow (quasi-static) invasion of a non- 
wetting phase. In terms of fuel cell operation, this simulation 
corresponds to the flow of liquid water (the non-wetting phase) 
from the catalyst layer through the GDL to the flow channel, via 
a path of the largest accessible pores. The algorithm for sim- 
ulating drainage in the network is as follows. First, an initial, 
low-capillary pressure is selected. The network is then scanned 
and all pore throats that could be penetrated at that given capil- 
lary pressure are marked as ‘open’, along with the pore bodies 
to which they are connected. Next, all distinct clusters of con- 
tiguous open throats and pores are found and labeled. Finally, 
all clusters that are connected to the injection face are identi- 
fied and are counted as invaded by the invading fluid. All pores 
and throats not connected to the injection face are returned to a 
‘closed’ state. In this way, the invading front of the non-wetting 
phase only reaches pores that are both topologically accessible 
from the injection face (i.e. through other invaded pores) and 
penetrable at the given capillary pressure. The algorithm pro- 
ceeds by increasing the capillary pressure in small increments 
and repeating the procedure until all pores and throats are open 
or filled with the invading fluid. The volume of non-wetting 
phase within pores that are invaded at each capillary pressure 
step is calculated and a capillary pressure curve is generated. In 
the present simulations, the injection of the non-wetting phase 
is always in the through-plane (z) direction. In terms of a GDL, 
the injection face is on one side of the paper and the exit face is 
the other side. 


2.6. Transport processes in the network 


2.6.1. Convection 

Determination of the flow rate and pressure drop across the 
pore network requires solution of the following mass conserva- 
tion equation over each pore: 


qi = >_gni(Pj — Pi) = 0 (5) 
J=1 


where i denotes the current pore, j denotes the neighboring pore, 
n is the number of neighbors, q; is the net flow through pore i, 
8h,ij is the hydraulic conductivity for flow between pore i and the 
neighboring pore j, while P; and P; are the pressures in each pore. 
The hydraulic conductivity, gn, of the pores and throats depend 
on their size and length and is determined from the following 


J.T. Gostick et al. / Journal of Power Sources 173 (2007) 277-290 283 


expression for square ducts [23]: 


2.28b4* 
2Lu 


gh = (6) 
where 2b is the size of the conduit opening, u is the fluid viscos- 
ity and L is the conduit length. L is equal to b for pore bodies and 
calculated for pore throats as discussed in Section 2.2. The total 
hydraulic conductivity for flow between two adjacent bodies is 
taken as the net conductivity for flow through half of pore i, the 
connecting throat and half of pore j. The hydraulic conductivity, 
gn, for each section is calculated using Eq. (6) and the net con- 
ductivity for the pore—throat—pore assembly, as shown in Fig. 3, 
is found from linear resistor theory for resistors in series: 
1 1 1 1 
= F + (7) 


&gh,ij &h,pi 8h,t &h,pj 


Eq. (5) is set up for each pore in the network to yield a system 
of linear equations that can be solved in conjunction with the 
prescribed boundary pressures on each side of the network to 
give the total flow (Q) across the network [17]. Once Q is known, 
the permeability of the network can be found from Darcy’s law: 


KA 
Q = — (Pin — Pout) (8) 
pl 


where K is the absolute permeability, Pin and Pout are arbitrarily 
chosen inlet and outlet boundary pressures. For flow in the z- 
direction, A = X YEA is the area of pore network normal to the 
direction of flow and / = ZLc is the length of the pore network in 
the direction of flow. X, Y and Z are the dimensions of the network 
in number of pores and Lc is the lattice constant, discussed in 
more detail in Section 3.2. 


2.6.2. Diffusion 

The diffusivity of the network is found in the same manner 
as for fluid flow. Fick’s law for binary diffusion of A through 
stagnant B is: 


cDap dxa cDap dxB d In xB 
= = cDAB 
l— xa dl xg dl dl 


where Dap is the binary diffusion coefficient, c is the mole con- 
centration, xa is the mole fraction of species A, xp is the mole 
fraction of species B (xg = 1 — xa), and / is the length of the 
domain. Using Eq. (9), the species conservation equation at each 
network node is then written: 


Na = 


(9) 


n 
ni = S gajn xp,j — ln xg) = 0 (10) 
j=l 


where n; is the mass transfer rate through pore i, xg j is the con- 
centration in the neighboring pore j, and xg,; is the concentration 
in pore i. gq is analogous to the hydraulic conductivity and is 
calculated for a given conduit as: 


cDap(2b)* 
ga = e (11) 


where Dap is the diffusion coefficient and 2b is the width of the 
conduit. The conductivity for diffusion through each half pore 


and throat is calculated using Eq. (11) and the net conductivity 
for the entire conduit is found from: 
1 1 1 1 
= +— + (12) 
Sdij Sdpi Sdt  Sdpj 


Upon solution of the system of species conservation equa- 
tions, the effective diffusivity of the network is found using 
Fick’s law: 


CD eA 
Nax eff 


(In XB,in — In XB,out) (13) 


where Deg is the effective diffusivity of the network. xB,in and 
Xpout are the inlet and outlet mole fractions of the stagnant 
species B. 


2.6.3. Multiphase transport 

In order to study conditions relevant to PEMFC operation, it is 
necessary to model the transport of gas and liquid as a function of 
water saturation in the GDL. This can be done by calculating the 
water and gas effective permeability and the gas diffusivity after 
the network has been partially invaded by the non-wetting phase 
(water), over a range of saturations. The general approach is to 
modify the conductivity of individual pores and throats as they 
become invaded by the non-wetting fluid and to recalculate the 
overall transport through the network. Since a certain amount of 
wetting phase is always present within pores and throats invaded 
by the non-wetting phase, due to late pore filling effects, careful 
attention must be paid to this modification, particularly in view 
of the fact that the precise geometry and connectivity of the 
remaining wetting phase is unknown. Two limiting cases are 
considered: 


Case 1. Once a pore is penetrated with the invading fluid 
(water), the residual wetting phase is no longer conductive. This 
case represents the most pessimistic scenario for gas transport 
since it leads to a highly obstructed and disconnected network 
with increasing invading fluid saturation. 


Case 2. The residual wetting phase within pores and throats 
invaded by the non-wetting phase maintains a connection with 
neighboring pores and offers limited conductivity to mass trans- 
fer through films and corners, which is modeled by assuming 
that the area for mass transport varies directly with the volume 
fraction of the conducting phase in a given pore. This case rep- 
resents the most optimistic scenario for gas transport since it 
neglects the tortuosity of the pore space containing the residual 
wetting phase. 


In general, for both cases the expressions for hydraulic and 
diffusive conductivity (Eqs. (6) and (11)) become: 


2.28b4 
Shi = Lu (Sop)” (14) 
and: 

cD(2b;)" 
gai = —— Sop)" (15) 


where Sop is the volume fraction of conducting phase in pore i. 
The exponents m and n control the behavior of the pore satura- 
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tion correction and depend on the conducting phase and case 
of interest. For Case 1, m=2 and n=1 for the non-wetting 
phase, while m and n are both equal to infinity for the wetting 
phase. The latter situation sets the conductivity to 0 for all pores 
that are invaded (swp < 1). For Case 2, m=2 and n= 1 for both 
phases. 


3. Model calibration 
3.1. Pore and throat size distribution 


The first step in the calibration of a pore network model is 
to identify the pore size distribution that enables the model to 
match experimentally determined drainage capillary pressure 
data. The computed drainage capillary pressure curves for SGL 
10BA and Toray 090 were compared to previously reported MIP 
data [12] for the displacement of air by mercury. Fig. 6 shows 
a comparison of the experimental data and the model curves 
obtained, while Fig. 7 shows histograms of pore size and throat 
size distributions used to generate these curves. The parameters 
for the Weibull distribution (Eq. (1)) obtained by fitting are listed 
in Table 3. The mean number averaged pore diameters for Toray 
090 and SGL 10BA obtained from these fit distributions are 19 
and 33 um, respectively. These values agree well with the results 
of Tomadakis and Robertson [24], who calculated pore size dis- 
tributions and mean pore sizes for solid models of various fiber 
arrangements and porosities. They also agree with similar data 
obtained recently by Schulz et al. [25] for simulated Toray 090 
and SGL 10BA materials. The fit in the high capillary pressure 
region obtained for the SGL 10BA sample was ignored since the 
pore space in this region represents sub-pore-scale roughness of 
the PTFE coating and binder materials (visible in Fig. 1b, ii). 


1.0 A MSP - Octane 
O MIP 
— —PNM - Octane 


—— PNM - Mercury 


o o o 
A O œ 


Non-Wetting Phase Filled Porosity 
oO 
iv 


10° 10° 10' 


relative frequency 


0.15 0.15 
0.1 0.1 
0.05 0.05 | 
0 0 il | 
0 10 20 30 0 10 20 30 


pore radius [microns] 


0.15 0.15 
0.1 0.1 
0.05 0.05 
0 0 

0 10 20 30 0 10 20 30 


throat radius [microns] 


0.15 0.15 
0.1 0.1 
0.05 0.05 
0 0 

(0 10 20 30 0 10 20 30 


throat length [microns] 


relative frequency 


relative frequency 


Fig. 7. Pore size, throat size and throat length histograms. (left) Toray 090 and 
(right) SGLIOBA. 
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Table 3 
Model parameters used for each material 
Toray 090 SGL 10BA 

Network size parameters 

Lc (wm) 25.2 40.5 
Pore size distribution parameters 

À 525 9 

K 3 3:5 

bmin (um) 5 9 

Xmax 0.95 0.9 
Late pore filling parameters 

s* 0.20 0.20 

n 1.00 1.00 
Throat constriction factors 

[ax, dy, œz] [1, 1, 0.9] [1, 0.95, 0.95] 
Pore correlation distances 

[Bx, By, Bz] [1, 1, 0] [2, 1, 0] 


The computed capillary pressure curves both rise more sharply 
than the experimental ones due to the use of a rather narrow pore 
size distribution, which is necessary to match the high porosity 
(see Section 3.2). 

To further assess the validity of the capillary pressure curves 
generated by the model, simulations were run with octane as the 
wetting fluid and air as the invading fluid. This corresponds to 
experiments performed using the method of standard porosime- 
try [12]. The advantage of considering this system is that octane 
is a highly wetting fluid and its contact angle can be confidently 
taken equal to 0°. It should be noted that the Weibull distribu- 
tion parameters listed in Table 3 and obtained above by fitting the 
model to the MIP data were also used for the octane—air system. 
The only parameters that differ were the surface tension and con- 
tact angle of octane (see Table 2). The good agreement between 
the simulated and experimental capillary pressure curves also 
shown in Fig. 6 supports the validity of the pore and throat size 
distributions selected. It is possible, however, that other pore and 
throat size distributions than those given in Table 3 could also 
lead to a match between the computed and measured capillary 
pressure curves. It is necessary to compare model predictions 
to other experimental results, such as absolute permeability 
and porosity, to improve confidence in the characterization of 
the two GDL materials in terms of the distributions given in 
Table 3. 


3.2. Lattice constant 


The lattice constant is the distance between pore centers in 
the cubic lattice. For a given set of pore sizes, adjusting the lat- 
tice constant controls the porosity of the network. For instance, if 
the lattice constant is large, then a significant amount of distance 
will exist between pores, thereby increasing the solid fraction 
and reducing the porosity. In the present work, the lattice con- 
stant was determined in the following manner. First, a pore size 
distribution was selected. Then an initial guess was made for the 
lattice constant and corresponding throat volumes (i.e. lengths) 
determined. This also allowed the porosity (£) of the network 


for a fixed total void volume to be calculated from: 
_ VHV 


e= (16) 
3 
fey og 


where Vp is the total pore volume of the network, V; is the total 
throat volume, X, Y and Z are the dimensions of the network 
expressed in terms of the number of pores and Lc is the lattice 
constant. The value of Lc was adjusted until the calculated poros- 
ity matched the experimental value for the material. Finally, it 
was verified that Lc was larger than the largest pore in the net- 
work to ensure that no pores overlapped. If this criterion was not 
met, then the pore size distribution was adjusted and the process 
repeated. 

Preventing the overlap of pores is necessary to avoid several 
inconsistencies in the network geometry, such as pore vol- 
umes being counted twice, throat lengths being negative and 
the center-to-center distance between pores being larger than 
Lc. Also, if pores were allowed to overlap, it would be trivial to 
match porosity, since any pore size distribution would suffice. 
Allowing such flexibility in the pore size distribution would also 
enable a near-perfect matching of the capillary pressure curve 
since an arbitrarily broad distribution could be used. On the 
contrary, requiring that no pores overlap tightly constrains the 
range of pore size distributions that can be used. For instance, 
if the pore size distribution is very wide, the network contains 
many small pores. Since the lattice constant is on the order of 
the largest pore, these small pores are surrounded by a substan- 
tial amount of solid, making it impossible to have a sufficiently 
high porosity. In the present work, it was necessary to use a pore 
size distribution that gave a slightly steeper capillary pressure 
curve than the experimental data (Fig. 6) in order to match the 
porosity. The ability to match the porosity, while still achieving 
a good agreement of the capillary pressure curves, is a strong 
indicator of the appropriateness of the pore size distributions for 
such high-porosity materials. 

The value of Lc obtained also indicates the appropriateness 
of the model geometry since Lc has units of length and repre- 
sents the spacing between pore centers. The lattice constant for 
Toray 090 has a value of 25.2 wm and indicates that 11 pores 
on average span the thickness of the material. SGL 10BA has a 
lattice constant of 40.5 um, corresponding to 10 pores across its 
thickness. These values are consistent with information on their 
structures obtained from SEM images of GDL cross-sections 
[12]. 


3.3. Absolute permeability 


The final aspect of the model calibration is to compare the 
permeability of the network with measured permeability values. 
Comparing the model results to permeability data allows for ver- 
ification of pore information that is not reflected in the capillary 
pressure curve, such as pore lengths and connectivity. It has been 
experimentally observed [21] that the in-plane permeability is 
higher than the through-plane permeability, a result that has been 
verified numerically [26] and analytically [27]. As discussed in 
Section 2.2.2, spatial correlation of pore sizes is included in the 
network in combination with slight throat constrictions in order 
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Table 4 
Transport results for each modeled material 
Permeability (x 10!? m?) Toray 090 SGL 10BA 

Experimental [21] Model Experimental [21] Model 
Ky 15 14 57 54 
Ky 15 14 45 48 
K; 9.0 9.5 37 39 
Effective diffusivity Numerical [28] Model Numerical [28] Model 
Dex 0.67 0.54 0.78 0.64 
Dey 0.67 0.54 0.75 0.61 
Dez 0.62 0.46 0.75 0.58 


to reproduce the observed anisotropy in the model. Measure- 
ments on Toray 090 indicate that the in-plane permeability is 
about 1.5-2 times higher than that in the through-plane direc- 
tion (Table 1). As discussed in Section 2.2.2, spatial correlation 
distances of [Bx, By, Bz] =[1, 1, 0] and throat constriction factors 
of [&x, @y, &z] =[1, 1, 0.9] have been used in order to fully match 
the permeability data. This procedure reproduces the anisotropy 
and gives good agreement between experimental data and model 
results, as can be seen in Table 4. The anisotropy of SGL 1OBA 
was somewhat more complicated due to the alignment of fibers, 
which caused the permeability to differ from one in-plane direc- 
tion to the other. To capture this, correlation distances of [B,, By, 
B-]=[2, 1, 0] are used along with throat constriction factors of 
[Æx æy, &z]=[1, 0.95, 0.95]. 


4. Model validation 
4.1. Effective diffusivity 


Determination of the effective diffusivity of the network pro- 
vides a useful means of independently verifying the chosen 
network geometry. Although experimental data for diffusion 
through GDLs are not yet available, limited numerical results 
have been presented by Tomadakis and Sotirchos [28] for fibrous 
materials with various arrangements of fiber alignment that cor- 
respond to GDL materials. The effective diffusivities predicted 
by the present model are compared with those of Tomadakis 
and Sotirchos [28] in Table 4. The agreement is reasonable con- 
sidering that no efforts were made to fit the model to those 
values. 


4.2. Liquid water injection 


Recent experiments have been performed by Benziger et al. 
[29] to measure the breakthrough pressure of liquid water in 
GDLs. In these experiments, the static pressure of a column of 
liquid water above a GDL is increased until liquid penetrates 
the sample. The pressure required for water breakthrough on 
various samples has been reported, including a sample of Toray 
120 with no PTFE treatment. This material is thicker than the 
Toray 090 considered here, but otherwise similar in structure. An 
experimental value of 3300 Pa was found, which compares with 
a value of 2483 Pa predicted by the present model. These values 


are within 25% of each other, which is reasonable considering 
that the materials are not necessarily identical. The reasonable 
agreement between the model and data suggests that the contact 
angle used for water on Toray 090 is reasonably correct. Similar 
data are not available for SGL 10BA. 


5. Results and discussion 
5.1. Relative permeability 


In the presence of two or more phases, the permeability of 
each phase P is reduced since the number of available pathways 
is reduced by the presence of the other phase(s). This effect is 
expressed in terms of relative permeability K;,p defined as the 
ratio of the effective phase permeability Kerr. p(sp) in the presence 
of another phase to the absolute permeability, or single phase, 
permeability K, i.e., 


Kepe.p(sp) = K Ky,p(sp) (17) 


where sp is the volume fraction of phase P in the network. Krp 
depends on the magnitude of saturation and history of saturation 
change (drainage or imbibition) and varies between 0 and 1. In 
studies employing continuum models the functional form of K, p 
has been assumed to be: 


K; p = Sp (18) 


where a is typically taken as 3 in the fuel cell modeling literature 
[30]. Eq. (18) is one of several empirical models of relative 
permeability and, to the best of our knowledge its applicability to 
two-phase flow in fuel cell materials lacks experimental support. 

Relative permeability calculations using the pore network 
model are based on the assumption that the pore-scale fluid occu- 
pancy is dictated exclusively by capillary forces — an assumption 
appropriate for low-capillary number displacements. To exam- 
ine the effect of GDL anisotropy, the effective permeability was 
calculated in the x, y and z directions through the network to 
yield the results plotted in Fig. 8. Non-wetting fluid invasion 
was always in the through-plane direction, which corresponds 
to liquid water flow from the catalyst layer through the GDL 
to the flow channels. Also shown in Fig. 8 for comparison are 
the curves obtained using Eq. (18) with a=3 for the two GDL 
materials. These results have been normalized for the intrinsic 
anisotropy of each material and so the directional differences 
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Fig. 8. Relative gas and liquid permeability as a function of water saturation in the network. (a) Toray 090 and (b) SGL 10BA. Both cases are shown for the gas 


relative permeability. Also shown is the result using Eq. (18) with exponent a=3. 


observed reflect the anisotropic effects caused by the presence 
of liquid water. This saturation-dependent anisotropy is due to 
the preferential spreading of the invading phase in the direc- 
tion of highest permeability, which is the direction of largest 
and most easily invaded pores. One of the major consequences 
of water spreading preferentially in the plane of the material 
is the significant reduction of gas transport in the through-plane 
direction. This suggests that the ideal GDL is one where the typ- 
ical anisotropy ratio is not only minimized, but reversed. Higher 
through-plane permeability would simultaneously limit detri- 
mental liquid water spreading and increase the intrinsic transport 
rates to the catalyst layer. A broad analysis of the effects of 
anisotropy in the GDL is given by Pharaoh et al. [31]. 

An important feature of these results is the non-zero liquid 
water saturation required for liquid water to break through the 
GDL. For Toray 090, the simulations show that liquid water 
saturations of 20% are necessary before a continuous liquid 
path spans the full thickness of the GDL. For SGL 10BA, the 
necessary liquid saturation is 10%. Below this critical liquid 
saturation, the liquid water permeability through the GDL is 
zero. This behavior is not described by the general form of the 
relative permeability function in Eq. (18) which predicts finite 
water permeability at vanishing water saturations. Nonetheless, 
the results obtained using Eq. (18) (i.e. the dashed line) are 
in rough agreement with pore network calculations of water 
relative permeability in the through-plane direction. 

Predictions of the relative gas phase permeabilities are also 
shown in Fig. 8. The gas phase permeability was calculated for 
both cases discussed in Section 2.6.3. In Case 1, the residual 
gas in an invaded pore offers no conductivity and gas flows 
entirely through the network of connected gas-filled pores. In 
Case 2, gas is allowed to flow through the non-filled portion 
of invaded pores. Both of these cases are somewhat unrealistic, 


for Case 1 prevents any flow through the space occupied by 
gas within water-invaded pores whereas Case 2 allocates to this 
space the hydraulic conductance of a straight conduit of reduced 
size. These cases, therefore, provide lower and upper bounds of 
gas permeability, respectively. The Case 1 results show that no 
gas conductivity exists above a critical water saturation of 65% 
for Toray 090 and 70% for SGL 10BA. A significant amount 
of gas still exists in the network at this critical saturation, but it 
is completely surrounded or trapped by the invading phase and 
is hydraulically disconnected from either the gas inlet or outlet 
face. Case 2 does not show a critical water saturation, since all 
trapped gas pores maintain some hydraulic conductivity. This 
case matches the behavior of Eq. (18) very closely. Since Case 2 
unrealistically allows gas transport to occur unimpeded through 
the corners of pores that are mostly filled with water, then Eq. 
(18) must be also be considered a limiting case. Eq. (18) requires 
a to be about 5 to match the model results for Case 1. 

Cases 1 and 2 exhibit other differences due to anisotropy 
as liquid water saturation is increased. Case 1 shows signifi- 
cantly reduced permeability in the through-plane direction due 
to spreading of liquid water in the x-y direction, whereas Case 2 
shows little to no anisotropy caused by additional liquid water. 
The latter effect arises because gas can leak through a pore even 
if it is mostly filled with water and allow pockets of trapped gas 
phase to contribute to mass transfer, thus minimizing the impact 
of in-plane liquid spreading. 


5.2. Dependence of effective diffusivity on water saturation 


The diffusion of gas from the flow channels to the catalyst 
layer is the predominant mode of reactant transport in a conven- 
tional PEMFC. As with gas convection, the presence of liquid 
water in the porous medium greatly reduces gas diffusivity. The 
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reduction of the diffusion coefficient due to the presence of liquid 
water is given as follows: 


Dest(é, Sw) = Dap Dyp(sw) f (£) 


where Degf(€,Sw) is the effective diffusion coefficient, Dap is 
the bulk diffusion coefficient, D,p(sw) is the relative diffusivity 
and f(¢) accounts for the reduction of diffusivity due to porosity 
and tortuosity. In fuel cell modeling literature, the Bruggeman 
approach is almost invariably adopted, leading to f(e)=e!°, 
although other estimates are available (see [32] and references 
therein), including one specifically for fibrous media [26]. Eq. 
(19) is analogous to Eq. (17). The function D,p(sw), which is 
here called relative effective diffusivity due to its analogy to the 
relative permeability, has not been as widely studied, particu- 
larly for GDL materials. Nam and Kaviany [11] have performed 
a numerical study using a rudimentary network model. The pore 
network studied by these authors lacked a pore size distribution 
and could not be tailored to specific GDL materials. More impor- 
tantly, in the model of Nam and Kaviany [11] water saturation 
was established with no regard for the physics of immiscible dis- 
placement. They suggested that the relative effective diffusivity 
decreases with the square of water saturation: 


(19) 


D,p(sp) = sp (20) 


where a= 2. Values of a= 1.5 are also commonly used [3] based 
on the assumption that the Bruggeman correlation for the effect 
of porosity also applies to the effect of liquid water saturation. 
The present model was used to calculate relative effective 
diffusivity in a GDL using invasion percolation concepts that 
more realistically simulate the configuration of water expected 
in an operating fuel cell. Specifically, liquid water was injected 
into the network in the through-plane direction to simulate liquid 
water flowing from the catalyst layer to the gas channels. The 
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present model also includes pore and throat size distributions 
that adequately reproduce both the absolute permeability and 
effective diffusivity through a dry network. The results are shown 
in Fig. 9 along with those using Eq. (20) with a=2. 

The difference between Case 1 and Case 2 is much more 
dramatic for gas diffusivity than for gas permeability. This is 
due to the fact that diffusional conductivity is proportional to 
the area available for transport, while hydraulic conductivity is 
proportional to the square of the area. Since the area for trans- 
port through a pore is assumed to be proportional to the volume 
fraction of a pore that is filled with gas, the diffusional con- 
ductivity is much less hindered by the partial filling of pores. 
The large discrepancy between these two limiting cases under- 
scores the need for experimental data concerning these transport 
processes. An argument against Case 2 is that not only does it 
fail to display a critical water saturation (above which effective 
gas diffusivity is zero), but it predicts significant diffusivity at 
near full-water saturation (Dra(sw =0.9)=0.1), which appears 
unrealistic. Case 1 shows a significant decrease in diffusivity as 
water invades the network. Compared to Eq. (20), diffusivities 
predicted by Case 1 can be several times lower. An exponent of 
a=5 would be necessary in Eq. (20) to approximate the behav- 
ior of the network model in this case. Clearly, current models 
could be significantly overestimating the transport rates through 
partially saturated GDLs. 

Also shown in Fig. 9 are the liquid phase diffusivities. These 
values are not of direct interest to PEMFC performance calcula- 
tions since liquid phase diffusion of reactants through the GDL 
is not significant. However, an area of research that is becoming 
increasingly active is the transport of ionic contaminants (e.g. 
Fe(II)), in the liquid phase. Thus, the presented results provide 
an estimate of diffusivities to be used in modeling contaminant 
transport in PEMFCs. 
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Fig. 9. Relative effective diffusivity as a function of water saturation in the network. (a) Toray 090 and (b) SGL 10BA. Both cases are shown for the relative effective 
diffusivity of the air phase. Also shown is the result using Eq. (20) with exponent a=2. 
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Fig. 10. Schematic diagram of modeled domain. x;,cy is the concentration of 
species i in the flow channel, x;,cy is the concentration of species i at the catalyst 
layer. 


5.3. Limiting current 


An effort was also made to use the present network model to 
predict the limiting current in an operating PEMFC assuming 
that the GDL is the sole source of mass transfer resistance. This 
was undertaken in order to determine if and when mass transfer 
resistance in the cathode GDL becomes a significant portion of 
the overall mass transfer resistance [33]. By estimating the max- 
imum rate of oxygen mass transfer that can be expected through 
a partially saturated GDL the limiting current was calculated 
and compared with typically observed values in operating cells. 

The modeled domain is shown in Fig. 10. The size of the 
domain is equivalent to 1 mm x 1 mm x ô, where ô is the GDL 
thickness (Table 1). This corresponds to a domain size of 
40 x 40 x 12 pores for Toray 090 and 26 x 26 x 10 pores for 
SGL 10BA. On the channel side of the domain, half of the inlet 
face is blocked to simulate the effect of 1 mm lands and channels. 
The conditions in the flow channel are taken as fully humidified 
air at 80°C and 10 kPa gauge. The catalyst layer is treated as a 
reactive interface where the oxygen concentration is zero (i.e. 
limiting current conditions). Since the cell is fully humidified 
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Fig. 11. Predicted limiting current densities as a function of GDL water satura- 
tion based on mass transfer through the cathode GDL. 


there is no water vapor diffusion and all water generated by the 
electrochemical reaction is in the liquid state. As a result, the 
mass flux through the GDL is considered to be molecular diffu- 
sion of Oz through a stagnant film of Nz and H20. This allows 
the multicomponent diffusion problem to be reduced to a binary 
diffusion problem, provided that the diffusion coefficient is cal- 
culated with appropriate consideration for the composition of 
the stagnant gas mixture [34]. Once the mass flux through the 
GDL is known, the current density is found from Faraday’s Law. 

The predicted limiting currents for both GDLs and both wet- 
ting phase conductivity cases are given in Fig. 11. The limiting 
currents through dry Toray 090 and dry SGL 10BA are very 
similar to each other. Although Toray 090 is 25% thinner than 
SGL 10BA, it is less porous and has a lower intrinsic effective 
diffusivity. These two factors offset each other and neither GDL 
is clearly better in terms of mass transfer performance under 
dry conditions. As water is added to the GDLs, however, the 
performance of the two materials diverges; the limiting current 
for SGL 10BA drops more quickly. This can be attributed to 
the increased spreading of liquid water in the x-y plane of this 
material. 

The overall behavior for both materials shows a dramatic 
decrease in limiting current as the GDL fills with water. At 
low water saturations (<10%), the predicted limiting current 
through the GDL is higher than in a typical fuel cell, which 
can be between 1 and 2 Acm7?. This indicates that at relatively 
dry conditions, the GDL is not the main source of concentration 
polarization, and performance is limited by other factors (i.e. the 
catalyst layer or electrolyte phase). When the GDL becomes wet, 
however, there is a significant reduction in the limiting current 
due to mass transfer resistance in the GDL. Case 1 predicts that 
at water saturations above 25% the maximum current density 
is less than 1 Acm~?, indicating that mass transfer resistance 
through the GDL could be a dominant factor limiting PEMFC 
performance. The limiting currents for Case 2 do not drop as 
sharply in the presence of water and 75% saturation must be 
reached before it reaches 1 Acm™?. 

At present, insufficient experimental evidence is available 
to fully understand the configuration and connectivity of the 
residual gas phase in GDL pores invaded by water. Some exper- 
imental evidence concerning the amount of liquid water in the 
GDL of an operating fuel cell does exist, however. Kramer et 
al. [35] used neutron imaging to measure the water content in 
the cathode GDL during fuel cell operation and found satura- 
tions between 25% and 35% at limiting currents between 0.6 and 
1.0 Acm~?, which corresponds very closely with the results of 
Case 1. Other neutron imaging studies suggest a limiting current 
above 1 A cm7? at somewhat higher water saturation (30-60%) 
[36,37], which lies between Case 1 and Case 2. Obviously, more 
conclusive evidence is needed to verify the present model, but the 
reasonable agreement with these experimental results does lend 
support to the applicability of the network modeling approach. 


6. Conclusions 


A pore network model was developed to help understand the 
multiphase flow properties of GDL materials and estimate their 
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multiphase flow and transport properties. A detailed description 
of the model was provided, with particular emphasis on integrat- 
ing into the model both qualitative and quantitative aspects of the 
microstructure of high-porosity fibrous GDLs. The model was 
calibrated to two commonly used GDL materials by adjusting 
the model parameters to match available experimental results, 
specifically the absolute permeability tensor and drainage cap- 
illary pressure curves. Material-specific relative gas and liquid 
permeabilities and diffusivities were computed as functions of 
water saturation under conditions of quasi-static drainage of air 
by water and transport rates through the pore network were deter- 
mined. Uncertainty regarding the configuration of the residual 
wetting phase (gas) in water-invaded pores of the material made 
it necessary to consider two limiting cases for gas transport: 
Case | in which residual gas phase is not conductive, and Case 
2 in which the conductivity of the pore space occupied by gas in 
water-invaded pores is optimal. The results of these simulations 
were compared with commonly used models of relative perme- 
ability and diffusivity. It was found that these models tended to 
agree with Case 2, which likely overestimates mass transfer in 
the gas phase. Alternative forms of these common models were 
proposed that match the pore network modeling results of Case 
1. This study further highlights an urgent need for experimental 
measurement of the effects of water saturation on water relative 
permeability and gas diffusivity. 

Limiting current calculations were performed by implement- 
ing PEMFC boundary conditions and physical parameters on the 
network model. The limiting current was estimated at various 
water saturation levels for a GDL section in which one-half was 
open to the gas channel and the other half was covered by a land. 
A dry GDL can support limiting currents of nearly 4 A cm7?, 
much more than is typically observed in operating fuel cells. 
When liquid water is present in the GDL, however, the predicted 
limiting current decreases rapidly to values typically observed 
in operating PEMFCs, indicating that mass transfer through the 
GDL may indeed be rate limiting at high current densities when 
the GDL is saturated with water. 
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